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We study dynamical mass generation in QED in (2 + 1) dimensions using Hamiltonian lattice 
methods. We use staggered fermions, and perform simulations with explicit dynamical fermions 
in the chiral limit. We demonstrate that a recently developed method to reduce the fermion sign 
problem can successfully be applied to this problem. Our results are in agreement with both the 
strong coupling expansion and with Euclidean lattice simulations. 



I. INTRODUCTION 



using large lattices. 
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Quantum electrodynamics in (2 + 1) dimensions 
(QED3) is a theory which shares a number of impor- 
tant features with quantum chromodynamics in 3 + 
1 dimensions (QCD) such as dynamical mass genera- 
tion I I g I I 0, |, |, and confinement 0. 
Since QED3 is super-renormahzable and has fewer de- 
grees of freedom than QCD, it serves as a valuable labo- 
ratory in which to test new methods and ideas related to 
these phenomena. Aside from its role as a testing ground 
for QCD, however, QED3 in itself plays an important 
role in solid state physics and in particular high-Tc su- 
perconductivity [1^. Recently several studies have 
pursued a new theoretical approach to cuprate supercon- 
ductors ^ in which one describes the phase tran- 
sition in the reverse direction, starting from the super- 
conducting state. In this picture the antiferromagnetic 
phase, for example, corresponds to spontaneous chiral 
symmetry breaking of massless two-flavor QED3. But 
there are also several other phases, and the large chiral 
manifold of degenerate states explains the complexity of 
the phase diagram. 

There are extensive studies using the Dyson-Schwinger 
equations |^ suggesting that chiral symmetry in QED3 
is dynamically broken if the number of fermion flavors is 
smaller than some critical number Nc ~ 3.3. However, 
the scale of this symmetry breaking (i.e. the magnitude 
of the chiral condensate) is extremely small, and there are 
also studies suggesting that chiral symmetry is bro- 
ken for all number of fermion flavors. Quenched lattice 
simulations have shown clear signs of chiral symmetry 
breaking, with a condensate {tpijj) 5 x 10~^ in units of 
e^, the dimensionful coupling constant Q. The situation 
for dynamical fermions however is not so clear, especially 
for an odd number of flavors. There have been Euclidean 
lattice studies in both compact and non-compact Q 
formalisms with different numbers of flavors, all suggest- 
ing a very small condensate. The most recent Euclidean 
lattice study of two-flavor non-compact QED3 suggests 
an upper bound for the condensate of '--^ 5 x 10~^ [lOl, 



On the other hand, Hamiltonian lattice studies of 
QED3 with one fermion flavor have suggested a rather 
large value for the chiral condensate. These studies were 
based on the strong coupling expansion |^ and varia- 
tional coupled cluster expansion 0. The obtained con- 
densate was ~ 0.28 1^, significantly larger than both 
quenched and two-flavor Euclidean lattice results and 
about two orders of magnitude larger than the Dyson- 
Schwinger results for one flavor QED3. 



In this paper we study chiral symmetry breaking in 
one-flavor massless QED3. To our knowledge our anal- 
ysis p7[ | represents the first non-perturbative simulation 
of lattice gauge theory in more than one spatial dimen- 
sion with explicit fermions. By explicit fermions, we 
mean that fermions are not integrated out to yield de- 
terminants of the Dirac operator. In the simulation pre- 
sented here, fermion dynamics are sampled explicitly us- 
ing fermion worldlines in a gauge-field dependent Hamil- 
tonian. From a theoretical point of view, this is an ideal 
framework in which to address the fermion structure of 
the ground state wavefunction. From a computational 
point of view, however, the approach presents profound 
difhculties such as the fermion sign problem and com- 
plex phase fluctuations due to the gauge field, both of 
which scale exponentially with the volume of the sys- 
tem. Therefore it is not likely that this approach would 
be possible for QCD in the near future. However, we 
do find that by employing the recently developed zone 
method 1 18 , we can control sign and phase problems suf- 
ficiently to study chiral symmetry breaking in massless 
QED3 on relatively small spatial lattices. 
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In our study, we find that in the strong coupling re- 
gion, y < 1, our results agree very well with the strong 
coupling expansion However, the agreement between 
the strong coupling expansion and our simulations breaks 
down around ?/ ~ 1, and for y > 1 we see a dramatic de- 
crease in the size of the condensate. These results are in 
agreement with Euclidean lattice simulations using stag- 
gered fermions ||^ , suggesting a very small condensate in 
the continuum limit, y 00. 
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II. QED IN (2 + 1) DIMENSIONS 

QED3 is a super-renormalizable theory, with a dimcn- 
sionful coupling: has dimensions of mass. This dimen- 
sionful parameter plays a role similar to Aqcd in QCD. 
In the chiral limit, it also sets the energy scale. We use 
4-component spinors, such that the fermion mass term 
is even under parity. With one massless fermion flavor, 
the Hamiltonian exhibits a global U{2) "chiral" symme- 
try. A fermion mass term breaks this symmetry to a 
C/(l) X U{1) symmetry. The question is: is this chiral 
symmetry broken dynamically? The order parameter for 
this symmetry breaking is the chiral condensate. 



A. Lattice Hamiltonian 

We start with the staggered fermion lattice Hamilto- 
nian on an Li x L2 spatial lattice B, 



physical = ^{We + Wb+Wf) 



with 



We 
Wb 
We 



r 

-M^(-ir+''^xt(r)x(rO 



(1) 

(2) 
(3) 



+ y^m{^)x\r)Uj{r)x{r + ]) +h..c. , (4) 



where ?7i(r) = (— ri2{r) = 1, y ~ 1/5^, and Up 
is the plaquette operator given by the product of Uj (r) 's 
circuiting the spatial plaquette anchored at r. 



Up{r) = Uiir)U2{r + i)Ul{f+ 2)t/T(f). 



(5) 



We use a dimensionless mass parameter /i = 2rn/e^ and 
a dimensionless coupling constant — e^a, where a is 
the lattice spacing. We also make the Hamiltonian and 
time dimensionless; the actual simulations are performed 
with the dimensionless Hamiltonian 



H 



2a 
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iJphysical = We+Wb+Wf, 



(6) 



in combination with the dimensionless time variable 

9' 



^ ■ 2^ ^physical; 

instead of the physical Hamiltonian, Eq. (|l|) 
We use the Dirac matrix representation 



io-1,2 
-icTi 2 



(7) 
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Assuming that Li and L2 are even, we stagger the 
fermion components at the four sites of a 2 x 2 unit cell 



XI (r) ~ (2ri,2r2) 

X2(f) ~ (2ri + l,2r2) 

X3(r) - (2ri + l,2r2 + l) 

X4ir) - (2ri,2r2 + l). 



(8) 
(9) 
(10) 

(11) 



In the continuum limit the staggered fermions correspond 
to one flavor of a 4-component fermion p. n9[. 
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For the states in our physical Hilbert space we choose 
a basis that is a tensor product of the gauge field and 
the fermion field degrees of freedom. For each gauge link 
field let us define the gauge field basis. 



(13) 



where each Xj is a real number in the interval [0, 27r). We 
let \X) be the tensor product of states |Xj(r)) at each 
link, 



\X) 



|X,.(rl) 



(14) 



Consider the Green's function 

Gx',a';X,a(At) - ( («' | ® (X' | ) e"^^* ( |X) ® |a) 

(15) 

where |X) (g) |a) and \X') \a') are two states in our 
physical Hilbert space, with \a) and la') general fermion 
states. If At is small and Xj{r) ~ ^'ji'r') for all r and j, 
then 



Gx' (At) 

« exp - 4St E (^^) - (^^)) ^-^tW, 

X (a'lexp [~AtW^] \a) , 



where 



Wb 
We 



(16) 

(17) 
(18) 



We can evaluate Gxf,a;;Xi,ai{t) for general initial and 
final states and arbitrary t by breaking the exponential 
in Eq. (|l^) into N equal time steps and inserting a com- 
plete set of states at each time step At = t/N . If At is 
small, the sum over intermediate states is dominated by 
consecutive states that are similar, i.e., the n"^ time step 



is dominated by states which satisfy X 



(n+l) 



X^"^ and 



J 3 

thus we can use Eq. ( |l6[) repeatedly. If we let Xi = 
Xf = X^^\ a, = and a/ = a^^\ then 



Gx/,Q/;Js:i,Qi(^) — J2 ^ 



X,a j^X 



(19) 



(1) 



where 



A 



(20) 



n=0 



X 
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X exp 



(21) 



In our simulation the gauge field configurations are 
updated using the Metropolis algorithm. For each new 
gauge configuration we compute the evolution of the cor- 
responding time-dependent Hamiltonian in the space of 
fermionic states. The sampling over fermion states is per- 
formed using the worldline formalism [ pO[ , which we now 
briefly discuss. 



B. Worldlines 

At the n^^ time step we have an exponential operator 
of the form 



where 



(22) 



Hr\r) 



H. 



(f) 



- "''-T'W r^xi^ (23) 
yJ^J^l^f}fll^x\r^x{r + i) + h.c.(24) 

yiV^.(rV4"'(^T ^t(^-^^(^-+ 2) + h.C.(25) 



In the following we use the shorthand "e" for even and 
"o" for odd values of r^, j — 1,2. Let us break up S^""^ 
into a product of four terms, 
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FIG. 1: Sample worldline configuration of a system with one 
spatial dimension. 



for a sample worldline configuration. For visual clarity 
the example we have drawn is a simpler system with only 
one spatial dimension. We have placed shaded squares 

(n) 

where the interactions S ■ , occur. In the case when two 
identical fermions enter the same shaded square we use 
the convention that the worldlines run parallel and do 
not cross. 

The sum over all worldline configurations is calculated 
with the help of the loop algorithm At each oc- 

cupied/unoccupied site, we place an upward/downward 
pointing arrow as shown in Fig. |^. Due to fermion num- 
ber conservation, the number of arrows pointing into a 
shaded square equals the number of arrows pointing out 
of the square. As a consequence of this conservation law, 
any valid worldline configuration can be generated from 
any other worldline configuration by flipping arrows that 
form closed loops. New Monte Carlo updates are there- 
fore implemented by picking a random closed loop and 
using the Metropolis condition to determine whether or 
not to flip the loop. 



C. Measuring the chiral condensate 



where 



S 



in) 



Each 'S'j"^/^ is the product of mutually commuting oper- 
ators which contain the interactions for an adjacent two- 
site system. With this decomposition of the time evolu- 
tion operator, one can trace out the worldline of any in- 
dividual fermion. In Fig. |l| we have drawn the worldlines 



exp 
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-E(>o(f)+Fr'(F) 



(26) 



.(27) 



For ^ = 0, the staggered lattice formulation reduces 
the chiral U{2) symmetry to a discrete symmetry gen- 
erated by a shift of one lattice spacing. To study chiral 
symmetry breaking on the lattice, we calculate the lat- 
tice condensate (V'V') i^i the chiral limit as a function of 
the lattice coupling y. The lattice condensate is related 
to the continuum condensate by the relation 



lattice 



1 



continuum 



(28) 



in the limit y oo. From here on (ipip) will denote 
We determine the lattice condensate by com- 
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D. Zone method 
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FIG. 2: Upward/downward arrows are drawn at each occu- 
pied/unoccupied site. 



puting the limit 



lim 



{{a\^{F\)e-"'{\F)®\a)) 



(29) 



where 



j_^(_l)n+..^t(^-.)^(^-.). (30) 



The state is a variational approximation to the gauge 
field ground state, 



\F) = J dXF{X)\X). 



(31) 



where 



F{X) = f,<^Y.r^os{Xi(r)+X2(r+i)-Xi(r+2)-X2(i^) ^ (33) 

and c is a real parameter we choose to optimize overlap 
with the true gauge field ground state . In our simula- 
tions we have used c = | , which appears to work well for 
both small and large y. The state |a) is the j/ = fermion 
ground state for /i > 0, a configuration where even sites 
are occupied and odd sites are unoccupied. The essential 
characteristic of the trial state \F) ® \a) is that it has 
non-zero overlap with the physical vacuum. 

In Eq. (|29|), for both numerator and denominator, the 
initial quantum state is the same as the final state. Any 
configuration of fermion worldlines can therefore be re- 
garded as a permutation of the initial fermions. Even per- 
mutations give a positive contribution while odd permu- 
tations come with a minus sign. Numerically, these mi- 
nus signs give rise to the fermion sign problem. With the 
worldline formalism we can keep track of these permuta- 
tions, and we use the recently developed zone method 
to manage the sign problem and well as phase oscillations 
due to the gauge field. 



The zone method |18 consists of introducing a special 
ni X n2 spatial sub-lattice or zone. See for example Fig. |^. 
We allow worldline configurations which may permute 



FIG. 3: Illustration of the zone method; an 6 x 8 spatial lattice 
with a 3 X 4 zone. The size of this zone is 17, as characterized 
by the number of links inside the sub-lattice. 

fermions lying inside this zone, but do not allow config- 
urations that permute any fermions lying outside of the 
zone. At intermediate time slices though the fermions 
are still allowed to wander through the entire lattice. In 
order to control phase oscillations associated with the 
gauge field, we use a different value of the coupling y' 
and/or fermion mass /i' when the fermions are outside 
the zone. We obtain physical results by extrapolation to 
the limit when the zone covers the entire Li x L2 lattice. 

As demonstrated in Ref. observables should scale 
linearly in the zone size provided that the zone size 
is larger than the characteristic "fermion wandering 
length" . For any finite values of the coupling y and of the 
time variable i, this wandering length is finite, even for 
massless fermions. Thus one can extrapolate the results 
from relatively small zones to the entire lattice. 



III. NUMERICAL RESULTS 

A. Zone Extrapolations 

There are different ways to define the size of a zone: by 
the number of lattice points or by the number of links in- 
side the zone. For large lattices it does not matter which 
is used. However, for the relatively small lattices we have 
used so far, it turns out that the best way to characterize 
the zone size is the number of links inside the zone. As 
an example, we show in Fig. ^ the lattice condensate for 
different zone sizes as function of the area of the zone, the 
number of points inside a zone, and the number of links 
inside a zone. A straight line fit to the condensate as 
function of the number of links gives a very good fit with 
a x^/d.o.f. of 0.6, whereas linear fits using the number of 
points or the area have a x^/d.o.f. of 6.1 and 8.5 respec- 
tively. Furthermore, the extrapolated result, using the 
number of links inside a zone, does indeed agree (within 
error bars) with the exact result. 
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FIG. 4: Numerical results for the condensate for different 
zone sizes within a 6 x 6 lattice, with fixed values of y — 0.4, 
y' = 0.1, n ^ fi' = 0, t = 1.5, and iV = 10. The total number 
of links in the lattice is 72; the total number or points is 36, 
as is the total area. The straight lines are linear fits to the 
data. 

To further test this method, we calculated the lattice 
condensate using dilFerent parameters y' and /i' outside 
the zone, while keeping the parameters y and /i inside 
the zone fixed. As can be seen from Fig. 0, the three 
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FIG. 5: Numerical results for different zone sizes within a 
6x6 lattice, for fixed values oi y = 0.5, fi — 0, t = 1.5, and 
= 10. The total number of links in the lattice is 72, where 
the three linear fits meet (within enlarged error bars). 

different sets of data points extrapolate to results within 
the error bars of the lattice condensate obtained without 
using the zone extrapolation method. The x^/d.o.f. of 
the linear fits are 0.3, 0.7 and 1.2 respectively, indicating 
that the numerical data are indeed on straight lines. 

Finally, in Fig. ^ we show results on a larger lattice. 
On an 8 X 8 lattice the method seems to work quite well, 



FIG. 6: Numerical results for different zone sizes within a 
8x8 lattice, for values of t/ = 0.4 and y = 1.0, both with 
fi — 0, t = 1.5, and — 10. The total number of links in the 
lattice is 128. 



although in this case we cannot compare our result with 
a simulation on the entire lattice. 



B. Finite size effects 

In order to avoid possible errors due to the zone ex- 
trapolation, we checked for finite size efi^ects without the 
zone method, which limits us to rather small values of y 
and coarse grids. In Fig. Qwe show for y — 0.4 the depen- 
dence of the chiral condensate on L, the spatial lattice 
size; on N, the number of time steps; and on t, the di- 
mensionless time variable. We see a slight dependence 
on L, which is actually smaller than our MC error bars. 
Note that the result for the 8x8 grid was obtained us- 
ing the zone extrapolation method, where the error bar is 
the error of the linear fit only. Within numerical error 
bars, our results are also independent of the number of 
time steps, N. 

The most significant finite-size effect is the dependence 
on t, as can be seen in the bottom panel of Fig. |^. An 
exponential fit of the type 

{^p^p){t) = = oo) + aoexp(-aii) (33) 

fits these data quite well for t > 1. However, for simula- 
tions at larger lattices and larger values of y the numerical 
errors are too large to do a proper finite-t extrapolation 
using such an exponential fit. 



C. Summary of main results 

In Figs. ^ and ^ we show our results for a range of val- 
ues of the coupling y. Most of the results are obtained 
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FIG. 7: Finite size eflFects for y = 0.4, ^ = 0: the condensate 
as function of the spatial square lattice of size LxL for A'' = 10 
and t — 1.5 (top), as function of the number of time steps A'' 
for L — 4 and t — 1.5 (middle), and as function of t for TV = 10 
and L — 4 (bottom). 



on 6 X 6 spatial lattices with = 10 using several differ- 
ent zone sizes. The error bars in Fig. || represent the 
error of the linear fit from our zone extrapolation. The 
error bars in Fig. ^ are our best estimate of the combined 
errors. They are dominated by the t — oo extrapolation, 
which is based on three (or more) different values of t 
where possible. For the largest values of y, y > 1.25, we 
could only establish upper limits for the condensate, due 
to the uncertainties in the t = oo extrapolation. 

Our results for the condensate indicate a dramatic 
change in behavior around y ~ 1: for y < 1 we agree 
within error bars with the 5/6 Fade approximant to the 
strong coupling expansion M . For y — > our results ap- 
proach the leading-order behavior in the strong coupling 
expansion 



y'^ii'ip) 



lattice 



0.5 y^ 



(34) 



as expected. However, for y > I we see a dramatic change 
in the behavior of the condensate, and a deviation from 
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FIG. 8: Our results for the lattice condensate as function of y 
for three different values of t on a 6 x 6 spatial grid, compared 
to Euchdean lattice simulations data from Ref. Q. 
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FIG. 9: Our results for the lattice condensate as function of 
y, obtained by extrapolating several different zone sizes for 
a finite value of t, compared to the strong couphng expan- 
sion [01. 



the strong coupling predictions: the value of the con- 
densate decreases rapidly. This strong decrease of the 
condensate with increasing y for 1 < j/ < 2 is in good 
agreement with the dynamical Euclidean Monte Carlo 
simulations by Burkitt and Irving j^] . 

It is not possible to determine at this time whether or 
not the condensate is small or exactly zero in the contin- 
uum limit, y oo. However, it is clear from our sim- 
ulations that the continuum condensate is significantly 



smaller than the prediction 



lattice 



0.284 



based on the strong coupling expansion. 
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IV. CONCLUDING REMARKS 

In this work we studied chiral symmetry breaking in 
one-flavor massless QED3, and our analysis represents 
the first non-perturbative simulation of lattice gauge the- 
ory in more than one spatial dimension with explicit 
fcrmions. While this approach is likely not practical for 
QCD in the near future, we were able to use the zone 
method to control sign and phase problems to study chi- 
ral symmetry breaking in massless QED3. 

We were able to resolve one puzzling issue regarding 
the size of the chiral condensate. Hamiltonian lattice 
studies had suggested a rather large value for the chi- 
ral condensate, whereas lattice simulations and Dyson- 
Schwinger studies indicated a value for the condensate 
about two orders of magnitude smaller. In our results we 
found that for y > I our results agree very well with the 
strong coupling expansion and the condensate appears 
to increase as y decreases. However for 1 < y < 2 we 
see a rather dramatic decrease in the condensate as y in- 
creases. These results are in agreement with Euclidean 
lattice simulations using staggered fermions 1^ . 



In future studies we would like to study the fermion 
structure of the ground state wavefunction and to com- 
pare and contrast what we see in the simulations with 
the coupled cluster variational state used in Ref . . We 
also plan to study the behavior of the chiral condensate 
as a function of fermion density. We note that stud- 
ies at finite density in this Hamiltonian formalism are 
no more difficult computationally than the simulations 
presented here. Since Euclidean lattice simulations of 
one-flavor QED3 at finite density are also afflicted by 
phase/sign oscillations (which make the computational 
effort scale exponentially with volume), comparison with 
explicit fermion simulations in the Hamiltonian frame- 
work could provide a valuable numerical check. 
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